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In field theory, particles are waves or excitations that propagate on the fundamental state. In 
experiments or cosmological models one typically wants to compute the out-of-equilibrium evolu¬ 
tion of a given initial distribution of such waves. Wave Turbulence deals with out-of-equilibrium 
ensembles of weakly nonlinear waves, and is therefore well-suited to address this problem. As an 
example, we consider the complex Klein-Gordon equation with a Mexican-hat potential. This sim¬ 
ple equation displays two kinds of excitations around the fundamental state: massive particles and 
massless Goldstone bosons. The former are waves with a nonzero frequency for vanishing wavenum¬ 
ber, whereas the latter obey an acoustic dispersion relation. Using wave turbulence theory, we 
derive wave kinetic equations that govern the coupled evolution of the spectra of massive and mass¬ 
less waves. We first consider the thermodynamic solutions to these equations and study the wave 
condensation transition, which is the classical equivalent of Bose-Einstein condensation. We then 
focus on nonlocal interactions in wavenumber space: we study the decay of an ensemble massive 
particles into massless ones. Under rather general conditions, these massless particles accumulate 
at low wavenumber. We study the dynamics of waves coexisting with such a strong condensate, and 
we compute rigorously a nonlocal Kolmogorov-Zakharov solution, where particles are transferred 
non-locally to the condensate, while energy cascades towards large wave numbers through local 
interactions. This nonlocal cascading state constitute the intermediate asymptotics between the 
initial distribution of waves and the thermodynamic state reached in the long-time limit. 
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I. INTRODUCTION 


Some nonlinear partial differential equations (PDEs) are simple enough and yet they describe deep and nontrivial 
universal processes in a broad range of physical applications. Well-known examples of such PDEs are the Korteweg- 
de-Vries and the Gross-Pitaevskii (GP) equations. It would not be an exaggeration to put into such a class a universal 
model that allows to understand the concept of spontaneous breaking of global symmetry. Consider for instance the 
following complex nonlinear Klein-Gordon equation, 

V’tt - A'*/) + (-1 + |V’|^)V’= 0, (1) 


where '0(x, t) G C. The linear part of this equation resembles the standard Klein-Gordon equation, where the square 
mass coefficient m? is replaced by — 1 . Such a negative square mass {w? = — 1 ) was originally believed to lead to 
superluminal particles. However, it was soon realised that such a field cannot sustain a localised particle-like excitation 
because the base state ip = 0 is unstable. 

Equation 0 conserves the total energy 


E = 


/ 


l^tp + ivV’p-iv.l^ 




dx, 


( 2 ) 


where the last two terms correspond to the famous Mexican hat potential U{'ip) = — + |'0|^/2; we therefore refer 

to equation Q as the Klein-Gordon Mexican Hat (KGMH) model. 

Equation is invariant to a uniform shift in phase of i.e., pj with a = const. The invariant corre¬ 

sponding to this 1 /( 1 ) symmetry is the “charge”. 


Q = \ J - Vdt'tp) dx, (3) 

where denotes the complex conjugate of pj. In this paper we focus mostly on the situation Q = 0, Q 7 ^ 0 being 
briefly described in section |mg 

In high energy physics, the KGMH model is a nonlinear a-model mm with a Mexican-hat potential. It is a simple 
relativistic model that describes the interaction between two fields: the first one corresponds to a massive particle, 
while the second one corresponds to a massless Goldstone boson. As such, the KGMH model provides a simplified 
description of a system of sigma mesons interacting with pions: the mass of the pions is neglected in this model, as 
it is much less than the one of the sigma meson. 

In this paper we study the KGMH system within the framework of wave turbulence, with direct numerical simula¬ 
tions (DNSs) to check the predictions. Wave turbulence O |4] deals with ensembles of weakly-interacting waves. It 
proceeds through the derivation of a kinetic equation that governs the time evolution of the spectral density of the 
wave field. Among its solutions are thermal equilibrium and generalized Rayleigh-Jeans distributions. However, the 
main use of the kinetic equation is to describe out-of-equilibrium situations. Indeed, for scale-invariant systems, one 
can compute out-of-equilibrium power-law solutions to the kinetic equation: these are the celebrated Kolmogorov- 
Zakharov spectra. More generally, the kinetic equation governs the evolution of an ensemble of waves (or particles), 
starting from an arbitrary weakly nonlinear initial condition. 

Such out-of-equilibrium initial conditions are encountered in particle physics and cosmology: they correspond 
for instance to heavy-ion collisions produced experimentally, or to the early universe. There has therefore been a 
significant recent interest in applying the wave turbulence approach in these domains lain]. 

Erom the point of view of wave turbulence theory, the KGMH model is a new object of study, which is interesting 
and important for several reasons: spontaneous symmetry breaking leads to two kinds of waves, or particles, with 
distinct dispersion relations. This contrasts with most wave turbulence systems, where one usually deals with a single 
kind of waves, like e.g. in the GP equation. In addition to transfers of energy from one scale to another, the KGMH 
system displays “reactions”: decay of massive particles into massless ones, and fusion of massless particles into massive 
ones. This brings new features to the wave turbulence dynamics, with strongly nonlocal interactions between different 
kinds of particles, that prohibit standard Kolmogorov-Zakharov cascades. 

Nevertheless, the KGMH model shares some similarities with the simpler GP equation: both have two positive 
invariants. Indeed, in addition to energy, the KGMH equation has an adiabatic invariant which is similar to the 
number of particles. Here, “adiabatic” means that this invariant is conserved exactly in the asymptotic limit of small 
nonlinearity only. This adiabatic invariance is sufficient to trigger wave condensation [7] , a phenomenon that we study 
both theoretically and numerically. In the particle physics context, this corresponds to Bose-Einstein condensation of 
pions ElIHlEl. 

The decay of the massive mode (cr-meson) into massless Goldstone bosons (pions) is the leading thread of this paper. 
In section]^ we propose a simple mechanical analogue of the KGMH model 0 , which provides physical insight into 
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FIG. 1. A lattice of magnetic pendulums: small magnets placed under the vertical position of each pendulum create the 
Mexican-hat potential. The interaction between neighboring pendulums results from dipole-dipole interaction and torsional 
springs between their fixation points. In the continuous limit this system follows the KGMH equation Q. 


the underlying dynamics. In section III we study the linear and weakly nonlinear dynamics around the minimum of 
the Mexican-hat potential (the uniform fundamental state). We compute the dispersion relations for the massive and 
massless modes, together with a reduced Hamiltonian which takes into account the dominant three-wave interaction 
between these modes. This dominant interaction involves two massless modes and one massive mode. We analyze 
it in section [IV| by deriving a simple set of ODEs for a single triad of interacting waves, which highlights the decay 
instability of a monochromatic massive wave into two massless monochromatic waves. 

In section [V| we develop the wave turbulence description of the KGMH model. We compute kinetic equations 
governing the evolution of the spectra of massive and massless modes. We identify two invariants of these equations: 
the energy, and an adiabatic invariant somewhat similar to a number of particles. The kinetic equations admit 
Rayleigh-Jeans equilibrium solutions, but no local Kolmogorov-Zakharov out-of-equilibrium solutions. We focus on 
such thermodynamic equilibria in section lYg For low enough energy per particle, the thermodynamic state displays 
wave condensation, the classical counterpart to the Bose-Einstein condensation (BEC) of pions. These findings are 
confirmed by DNSs, which show unambiguously that BEC is observed over a very long time, even in systems where 
the number of particles is only an adiabatic invariant. 

The out-of-equilibrium dynamics involve strongly nonlocal processes which we consider in sections m and |vmi 
The first section addresses the decay instability for isotropic ensembles of massive waves: on the one hand, the kinetic 
equations allow to predict the distribution of the massless waves that appear during the linear stage of the instability. 
On the other hand, the thermodynamics of section |Vl| predict the the reaction yield of the decay instability: starting 
from massive-waves (a-mesons) at a single wavenumber, we predict the fraction of these waves that has decayed into 
massless waves (pions) in the long-time limit. 

Section |VHI| finally addresses the behavior of ensembles of waves coexisting with a strong condensate of massless 
particles. We derive rigorously a reduced kinetic equation that describes a nonlocal transfer of particles to the 
condensate, together with a local Kolmogorov-Zakharov energy cascade towards high wave numbers. We see this 
cascading state as an intermediate asymptotics that describes the dynamical process of condensation: it makes the 
link between the initial distribution of particles and the final thermodynamic condensed state. 


II. AN ANALOGOUS MECHANICAL SYSTEM 

One can gain some insight into the physics described by equation 0 by considering a system from classical mechanics 
that follows the same equation. An example of such system is the square lattice of pendulums sketched in figure 
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The pendulums can move in two angular directions, ( and The corresponding kinetic energy of a pendulum is 
J(C^ + ^^)/2, with J the moment of inertia of a pendulum. Each pendulum consists of a magnetized bead attached 
to a rigid rod, the magnetic dipole of the bead being alined with the rod. Neighboring pendulums interact through 
two mechanisms: first, two dipoles repulse each other. Second, we assume that the fixation points of the pendulums 
are connected by torsional springs. Dipole-dipole repulsion efficiently propagates longitudinal perturbations through 
the lattice, while torsional springs propagate transverse oscillations. The stiffness C of the torsional springs is chosen 
to match the strength of the dipole-dipole interaction, so that the interaction energy between the pendulums reads 

f E (o-Ci)^ + te-ei)^ (4) 

(i,ji) neighbors 


where and are the angular degrees of freedom of the n-th pendulum, and we consider dipole-dipole interactions 
between nearest neighbors only. We further assume that the relative displacements of the pendulums are small and 
consider the lowest order description of the dipole-dipole interaction. 

Finally, to create the Mexican-hat potential we place a magnet under the rest position of each bead. Gravity drives 
a given pendulum towards the vertical, but the repulsion between the bead and the magnet displaces the equilibrium 
position away from the vertical, at a distance r = vq from the vertical axis. The combined effects of gravity and 
of the magnet result in a Mexican-hat-shaped potential V{r). If the length i of the rods is much longer than ro, 
we can use the approximate relation ^ f^e neighborhood of ro, we approximate the potential by 

V{r) ^ —V 2 r^ + The Lagrangian of the system follows from the subtraction 

of the Mexican-hat potential energy and neighboring interaction term to the kinetic term. 


^ — 2 + Cn ) - V4i^{Cl + -h E2^^(Cn + Cn) 


f E - 0)' + • 

(ij) neighbors 


( 5 ) 


We consider the continuous limit and introduce the complex variable '0(x, y, t) = C(x, t). After a rescaling 

of time, space, and this Lagrangian takes the following dimensionless form. 


£ = y - iV'i/jp + (6) 

which corresponds to equation ([T) in two spatial dimensions. 

In this analogous system, the charge Q corresponds to the sum of the vertical angular momenta of each pendulum 
around their fixation points. A nonzero value of Q means that the initial condition (and the subsequent evolution) 
favors one direction of azimuthal rotation of the pendulums. 


III. WEAKLY NONLINEAR DYNAMICS AROUND THE MINIMUM OF THE POTENTIAL. 

A key feature of the KGMH model is that it displays spontaneous U{1) symmetry-breaking: the fundamental state, 
defined as the state of lowest energy, is not the usual ip = 0. Instead, it is a uniform state ip = ipo = const., with 
|'0o| = 1: Ip is uniform, time-independent, and lies in the minimum of the Mexican-hat potential. We wish to study 
the dynamics of the system around this fundamental state. 

Up to a phase-shift in the definition of ip, we can choose ipo = 1 and consider perturbations of this state, 

'0(x,t) = l + 0(x,t). 

Substitution into equation Q yields 

(ptt ~ ^4^ -|- 0 + 0* + 200* -|- 0^ + I0p0 ~ 0 • C^) 

We focus on the weakly-nonlinear dynamics arising from the quadratic terms of this equation, and therefore neglect 
the cubic nonlinearity. In variables A = 5R(0) and y = ^(0) we have 

\tt — AA + 2A + 3A^ T = 0? (8) 

Xu - Ay + 2Ay = 0. (9) 

The linear parts of these two equations provide two dispersion relations. For a plane wave proportional to 
they give respectively 


wg = V 2 + fc 2 , 
Wk = k, 


( 10 ) 

( 11 ) 






5 


3 



FIG. 2. Two types of waves propagate in this system: branch a corresponds to radial oscillations of the pendulums, whereas 
branch b corresponds to azimuthal oscillations of the pendulums, along the minimum of the Mexican-hat potential. 


where k = |k|. We refer to these two dispersion relations as branch a and branch b. They are plotted in figure 
Branch a has a nonzero frequency for vanishing wavenumber, 7 ^ 0 for /c = 0. In terms of particle physics, this 
mode has a nonzero rest energy and thus corresponds to massive particles. By contrast, the frequency of 6 -waves 
vanishes for k = 0: branch b corresponds to a massless Goldstone boson. This is a basic model for coexisting sigma 
meson (branch a) and pion (branch b) fields, in which the pion’s mass is neglected. 

In the lattice of pendulums sketched in figure the massive mode corresponds to oscillations of the pendulums 
along their radial direction. In this mode, each pendulum feels two restoring forces, one from the potential well V (r) 
and one from the coupling to its neighbors. By contrast, mode b corresponds to oscillations of the pendulums along 
their azimuthal directions. The restoring force originates only from the coupling with the neighboring pendulums, 
and if all the pendulums move in phase along their azimuthal direction (6 mode with infinite wavelength) then there 
is no restoring force: mode b is massless. 


A. Hamiltonian formulation 


The system (| 8 ][^ follows from the following Hamiltonian, 


\p^ + + ^(VA)2 + A2 + I(Vx)' + A3 + Ax' 


dx, 


( 12 ) 


where p and q are the momenta conjugate to the variables A and x respectively. Hamilton equations are 


At — r )-Pt — 
op 

5H 

Xt yQt 

Sq 


“a ’ 

Sx ' 


(13) 


( 14 ) 
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We consider these equations in a d-dimensional periodic cube of side L and introduce Fourier series, 

/ \ 

X(x,i) 

p(x,t) 

V q{^,t) ) 

where the sums are over k G ^Z'^. Because these four fields are real-valued, (A_k, X-k,P_k, ^-k) = (^k> Xk;Pk;‘?k)- 
In Fourier space, Hamilton equations therefore become 


= E 


( Ak(t) \ 
Xk{t) 
Pk{t) 

V 9 k (f) J 


^ikx 


(15) 



SH SH 

(16) 


, SH SH 

Xk - ^,9k - 

(17) 

with 


H — H 2 + Hint^ 

(18) 

where 


H2 = \Y. [iPkl" + (^£)^|Ak|^ + 141^ + (c^^)^IXkl^' , 

k 

(19) 

and 


Hint — ^ ^ l^-^ki'^k2'^k3 H“ '^kiXk2Xk3 ^(^1 + \^2 H“ ^ 3 ), 

(20) 

where 5 is 

ki,k2,k3 

a Kroenecker symbol, i.e., b(k = 0) = 1 and b(k / 0) =0. 



B. Normal variables 


Let us introduce the normal variables Uk and 6k that diagonalize H 2 , 

C + <^Lk - V^^(^k — ^-k) 

^k = -^=^,Pk = - 


Xk = 


6 k + 6 ^j^ 


,^k = 


iV2 

- &*-k) 

iV2 


We refer to Uk and 6k as the respective amplitudes of a- and 6-waves (or modes, or particles). 
In these variables the equations of motion are 

5H 


.. SH 


where H = H 2 ^ Hint, with 


if2 = EKl«k|A4|6k( 


( 21 ) 

( 22 ) 


(23) 

(24) 

(25) 


and 


ki,k 2 ,k 3 


1 


A ,a 


— E 


rCo’ 2 Co’^ 

1 




= [(aia2a3<^(ki -f k2 -f ka) -f 3aia2a36(ki -f k2 - ka) 
[ai626a<5(ki -fk2 -l-ka) 


2 \/ 2 . . . 

ki,k2,k3 

-|-2ai62636(ki -f ka - ka) -f a\b2bz8{ki - ka - ka)] -f c.c.. 


( 26 ) 
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where c.c. stands for “complex conjugate”, and we introduced the shorthand notations bi = b]^. 

and cj? = ujt.‘ 




The terms in the Hamiltonian (26) correspond respectively to 3-wave processes of type a + a-l-a ^ 0,a-l-a ^ 
a, a + 6 + 6 ^ 0,a + 6 —> b and 6 -f 6 ^ a. In a weakly nonlinear description of the system, the only relevant 
nonlinearities are those for which both a wavenumber and a frequency resonance conditions can be satisfied. The 
other nonlinearities produce non-resonant terms only, which do not contribute to the weakly nonlinear dynamics. 
We show in appendix [A| that, out of the processes outlined above, the only one allowed by the frequency resonance 
condition is b ^ b ^ a. We can therefore discard all the other terms in as they do not contribute to the weakly 

nonlinear dynamics, and simply write 


Hint= y] ^23«1^2&3<^i*2 + C-C., (27) 

ki,k2,k3 

where we have switched to a shorthand notation of Kroenecker deltas, — k2 — ks), and introduced the 

interaction coefficient 




1 1 

2\/2 


The equations of motion become 


dk = -iWkflk - * bib25i2, 

ki,k 2 

k = y] Vkaib^dk. 

ki,k 2 


(28) 


(29) 

(30) 


A similar Hamiltonian system with two types of interacting waves was studied in m- However, in that study the 
dominant process was of the type a -\- b ^ 6, which leads to very different dynamics than the process b b ^ a 
considered here. 


C. Finite-charge Q 7 ^ 0 


We have considered so far only perturbations around the minimum energy state, which is uniform ip lying in the 
minimum of the Mexican-hat potential. This corresponds to a vanishing charge invariant, Q = 0. We now briefly 
consider nonzero values of Q. Writing the field in polar coordinates, pj = r(x, the charge reads 


Q = 


/■ 


^0 dx. 


(31) 


A nonzero charge therefore indicates a preferred direction of rotation in the complex plane. Uniform states with 
nonzero charge can be sought in the form '0o(^) = where R and ujq are real constants. Substitution into the 

KGMH equation ^ leads to 


R^ = 1 


^0 5 


(32) 


i.e. there is a continuous family of uniform solutions to the KGMH equation, corresponding to different values of the 
charge invariant Q = R^ujqL^. 

Once again, the lattice of magnetic pendulums provides a simple interpretation for these states. From equation 


(31) the charge is the sum of the vertical angular momenta of the pendulums about their fixation points. The uniform 


solution with cjo = 0 corresponds to all the pendulums parallel and at rest in the minimum of the Mexican-hat 
potential, whereas uniform solutions with cjq 7 ^ 0 correspond to all the pendulums parallel and rotating at the same 
angular frequency ujq around the vertical. For the latter solutions = 1 + cJq > 1- the field is not exactly in the 
minimum of the Mexican-hat potential, but instead it is “bobsleighing” on its outer edge. In the system of pendulums, 
this “bobsleighing” results from the centrifugal force acting on each pendulum. 

To study perturbations around the uniform state with nonzero Q, we write pj = '0o(t)(l-F(/)(x, t)) = 1 /I -h + 

0(x, t)), and consider infinitesimal perturbations 0 <C 1. Substitution into the KGMH equation gives at linear order 
in 6 


<Ptt + ~ A0 -h (1 + ^o)(^ ^*) ~ ^ 


(33) 







and after taking the real and imaginary parts, with 0 = A + iy, 


(34) 

(35) 


Xu — AA + 2(1 + cjq)A — 2cJoXt — 0 ? 

Xu - Ax + 2uJoXt = 0 . 


As compared to the situation (5 = 0 , the fields A and x ^-re now coupled at linear order by terms proportional to ujq. 
Inserting a plane wave structure for these two fields and asking for non-trivial solutions yields the following dispersion 
relation 


— 2 uj ^ { k ^ -|- 1 “I" 3cjq) a { k ^ -1- 2 -|- 2cjq) — 0 , (36) 

which has two branches of solutions, 

1 + ± 5(1 + 3wg)^ + , (37) 


In the limit cjo ^ 0, we see that branch corresponds to a-modes, whereas uj- corresponds to b-modes. For cjo 7 ^ 0, 
the mode remains massive, with lim/c^o <^+ 7 ^ O 5 and the uj- branch remains massless, with lim/c^o = 0- Note 
that the speed of the \ow-k 6 -particles is now reduced, duj-fdk = y ^(1 + cJq)/(I + 3co’o)? and the a-particle’s mass is 
+ 3 co’q times greater than in the cjq = 0 case. At the linear level, the structure of the problem is nevertheless quite 
similar to the Q = 0 situation. 

Interesting new phenomena appear at the nonlinear level. As for the Q = 0 situation, one can compute normal 
variables Uk and 6 k diagonalizing the quadratic part of the Hamiltonian, Uk (resp. 6 k) corresponding to Cc;+ (resp. 
UJ-). These normal variables are now superpositions of the A and x fields, together with the corresponding momenta: 


Ak 


Ak 

Pk 

& 


ttk + cj?. — k‘^ 6k — 61j^ 

ak-a*_k 1 6k + 6^k 
^ i72 \/2 


- k"^) ak - a;!^ 

71 ^ 



ttk + 

~W~ 


+ \/6jr 


6k + 61t, 

^2cj_cJq 2 
^k-6*_k 
^^/2 ’ 


(38) 

(39) 

(40) 

(41) 


The weakly nonlinear limit is still governed by three-wave interactions, but the dispersion relations (37) now allow 
for new processes: a + a ^ a and a + 6 ^ 6 are still impossible, but, in addition to 6 + 6 ^ a, the processes a ^ a + 6 
and a + 6 ^ a are compatible with the dispersion relation. This possibly triggers interesting new dynamics. However, 
because the normal variables Uk and 6 k have rather intricate expressions when (5 7 ^ 0 , the computation of the kinetic 
equation is more involved than for (5 = 0. The wave turbulence analysis of the present paper therefore focuses on 
(5 = 0 , and we leave the case (5 7 ^ 0 for future work. 


IV. DECAY OF THE MASSIVE MODE INTO TWO GOLDSTONE WAVES 

The dominant three-wave process of the weakly-nonlinear KGMH dynamics can be highlighted by studying a single 
triad of interacting modes. Consider the dynamical equations in the weak-nonlinearity limit, ^ and (§, and expand 
the fields as 


A = eAo(t,T) + e2Ai(i,T) + 0(e=*), (42) 

X = exoit,T)+e\i{t,T) + Oie^), (43) 

(44) 

where the slow time scale is T = et. To order e, equations ^ and ^ reduce to 

dtt^o ~ AAo + 2Ao = 0, (45) 

dttXo - Axo = 0. (46) 
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In Appendix we saw that triads of interacting mode are made of two b modes and one a mode. We consider a 


single such triad in the solution to the order 0{e) equations: 

Ao = + c.c., (47) 

Xo = + c.c.. (48) 

To obtain non-trivial dynamics at next order, we further assume that the resonance conditions \/fc^ + 2 = ki + k 2 
and k = ki + k 2 are met. To order e^, the equations become 

dtt^i ~ ^Ai + 2Ai = —2dtT^o — SAq — Xo ) (49) 

dttXi - Axi = -25 (tXo - 2AoXo • (50) 


The solvability condition requires the right-hand side of these equations to have no terms that are resonant with the 
operator on the left-hand side. We insert solutions (47) and (48) into equations (Eojand (50), before demanding that 

the rhs of (49) has no term proportional to of^5Q| has no terms proportional to 

g 2 (/eit-ki x) g 2 (/c 2 t-k 2 x)^ These three constraints lead to the following system of ODEs for the slow evolution of 

the wave amplitudes 


(ItA = 


2^/W 


= ^ 1^2 : 


drB, = , 

Ki 


= -^ABl . 

A 2 


(51) 

(52) 

(53) 


Consider a single 6 -wave: {A = B 2 = 0 , Bi 7 ^ 0 ) is a time-independent solution to the system of ODEs. We study 
the stability of this 6 -wave by considering infinitesimal perturbations A ^ 1 and ^2 ^ 1 - The linearized evolution 
equation for these perturbations is 


(IttF + 


\Bi? 

2k2\/WT2 


F = 0, 


(54) 


where F = A or F = B 2 . The b-w&ve is stable, in the sense that the perturbations of A and B 2 oscillate but remain 
small. 

By contrast, if one starts off with a single a-wave (A 7 ^ 0 , = ^2 = 0 ) and perturbs it with infinitesimal Bi and 

B 2 , the subsequent evolution of these perturbations is governed by 


(IttB 1-2 


kik: 


Bi.2 — 0 , 


(55) 


which admits some exponentially growing solutions: the a-wave is unstable and decays spontaneously into two 6 -waves. 
In terms of particles physics, this instability would describe the disintegration of a massive particle (e.g. the cr-meson) 
into two massless Goldstone bosons (two pions). 

The growth rate 7 = \A\/y/kik 2 of this instability is a first indication of n onlo cal interactions in Eourier space: 
recalling the resonance conditions, 7 is maximal when kik 2 = 1/2 (see condition (A 2 )), i.e., when ki = (\/2 + k‘^-\-k)l2 
and k 2 = (V2 + k?‘ — k)/2 (or vice-versa). Nonlocal interactions occur in at least two cases; if A: <C 1, then ki ~ 
A :2 V2 » fc. If fc ^ 1 , then either ki k or k 2 ^ k. Such nonlocal interactions are an important feature of wave 
turbulence in the KGMH model. 

To study the subsequent evolution of this instability, we solve equation ^ numerically. We focus on 2D numerical 
simulations, which corresponds to the dimension of the equivalent lattice of pendulums. We truncate the system at 
large wavenumbers. Such a truncation naturally arises from the discreteness of the lattice of pendulums. Eor more 
general applications of equation ([^, the truncation can be seen as a crude modelization of the ultraviolet cut-off 
arising from Bose statistics. 

Note that the domain size is an important parameter of this problem: the dispersion relation is not self-similar, 
with non-relativistic waves for /c ^ 1 and relativistic ones for k ^ 1. To accurately describe a mixture of relativistic 
and non-relativistic waves, we need to choose a domain size that is much larger than 27r. 

The numerical code is described in appendix|^ it uses a standard pseudo-spectral method, with an exact integration 
of the stiff linear wave operator and Adam-Bashforth time-stepping to deal with the nonlinear terms. Dealiasing is 
performed using the one-half rule for this equation with cubic nonlinearity. 
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FIG. 3. Destabilization of a small-scale massive wave into massless ones. The left-hand panels are the A helds, the middle 
panels are the x fields, and the right-hand panels are the ID spectra of A (blue) and x (red-dashed) as a function of k. Here 
L = 207r, and from top to bottom: t = 0, t = 30, t = 60, and t = 120. The unstable triads involve a small-scale massless mode 
and a large-scale one. 


In this section the domain size is [0,207r]^ with a resolution of 256^. After dealiasing we keep 128 modes in each 
direction: k G x {—64; 64}^. In figure|^we show snapshots of the A and x fields, together with their one-dimensional 
spectra, i.e., their spectra integrated over the direction of k. The initial condition is a small-scale progressive a-wave 
with k = (0;3), plus some weak background noise. The instability rapidly sets in and 6-waves appear. As can be 
seen in equation ( [^ , the maximum growth rate of the instability is achieved when the product kik 2 of the two b 
wavenumbers is small. Minimizing this quantity while satisfying the resonance conditions suggest that ki 1 and 
k 2 ^ k. Indeed, in figure we observe the emergence a combination of large-scale 6-waves, together with 6-waves 
around the wavenumber k of the initial a-wave. 

By contrast, in figure we show similar snapshots for a simulation initiated with a large-scale a-wave, with k = 
(0; 0.3). Once again, 6-waves rapidly appear as a result of the decay instability. However, for such a large-scale a-wave, 
the frequency resonance condition imposes that one of the 6-waves has a frequency (and therefore a wavenumber) of 
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FIG. 4. Destabilization of a large-scale massive wave into massless ones. The left-hand panels are the A fields, the middle 
panels are the x fields, and the right-hand panels are the ID spectra of A (blue) and x (red-dashed) as a function of k. Massless 
waves appear at the crossover scale /c ~ 1. Here L = 207r, and from top to bottom: t = 0, t=120, t = 270, and t = 570. 


order unity. But then the resonance of the three wave vectors imposes that the second 6-wave also has a wave number 
of order unity and almost antiparallel to that of the first 5-wave. As a consequence, in figure the instability sets in 
with 5-waves having wave numbers of order unity. 

In the following we extend the description of these nonlocal interactions to random ensembles of waves, using the 
framework of wave turbulence. 


V. WAVE TURBULENCE 


We now consider ensembles of weakly-nonlinear waves with homogeneous statistics in an infinite domain. For such 
systems wave turbulence can be used to derive wave kinetic equations that govern the evolution of the wave spectra. 
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A. Wave kinetic equations 


Let us define the wave action spectra of the modes a and h\ 


L 


= (^) (KP). 


= (^) 


here the angular brackets denote an ensemble average over many realizations. 

In appendix we apply standard wave turbulence techniques to equations (29) and (30), which leads to the 
following kinetic equations for the evolution of the two spectra, 




47r J IV 2 P (nlnl - dk,dk2, 

Stt j |y2kP««2 + <«k-'«-2«k)<^(‘^12k)<^Lf^kidk2. 


(56) 

(57) 


B. Conservation laws for the kinetic equations 


Kinetic equations (56) and (57) conserve the wave energy, 

£ = L'^ j (wgnk + w^n^) dk, 


as well as the following particle invariant, 


J\f = L' 


'/(2< + 4) 


dk. 


(58) 


(59) 


This second invariant is nontrivial, and it can be guessed only once the dominant three-wave process is identified. 
Indeed, the process a ^ 6 + 6 is similar to a chemical reaction where a is a “molecule” made of two b “atoms”, and 
J\f would be the total number of such b “atoms”: 1 per 6-wave, and 2 per a-wave. JV is also the number of b particles 
that would be in the system if the reaction a ^ b b were complete. For brevity, in the following we simply refer to 
Af as the number of particles. 

An important comment to make here is that Af is only an adiabatic invariant: it is conserved by the weakly 
nonlinear dynamics, but there is no similar invariant for the strongly nonlinear KGMH equation. By contrast, £ is 
the quadratic part of the global invariant the latter being conserved even in the strongly nonlinear regime. 


C. Thermodynamic equilibrium 


One can easily check that the kinetic equations admit solutions corresponding to energy equipartition in k-space, 

T 


= 






(60) 


for any temperature T = const m, and solutions corresponding to an equipartition of the particle invariant, 

= 2n^ = const. (61) 

More generally, the equilibrium Rayleigh-Jeans (or thermodynamic) solution is given by 

T T 

(62) 


+ 2/i 




where = const is a chemical potential. 

Finally, we remark that any initial condition with a vanishing spectrum = 0 and an arbitrary spectrum is 
a stationary solution to the kinetic equations (56) and (57). One can also show that spectra with 6-waves at low 


wavenumber only are also steady solutions to the kinetic equations. Such distributions with only a single type of 


particles are described in details in sections VII and VIII 
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FIG. 5. Integration domain Va{k) for k < 1/V^ (left-hand panel) and k > 1/V^ (right-hand panel). 


D. Isotropic systems 


Consider the class of isotropic spectra, = 
variables describing the directions of ki 


and = 


k = |kL and integrate out the angular 
and k 2 in the kinetic equations (56) and (57). Because ^ and 


where 


6{uj 


^ 2 ) ^ 


abb \ 
kl2) 


-yC 


variables is (5^2- This leads to 


are independent of the directions of ki and k 2 , all that needs to be integrated over the angular 


II 

f a 2n\nl) (kik 2 Y ^ dkidk 2 , 

(63) 


h^(k) 


00 

II 

f a i n\n\) {kik 2 Y ^ dkidk 2 , 

(64) 


h,(k) 



where 


and 


Sl2k — 


1 

Z/'W Zjyrv rv^ 1 rv 1 rv^rvQ^J rv rv^ 


for d = 2 , 


(65) 


Sl2k = 


1 


2kkik2 


for d = 3 . 


( 66 ) 


The integration domains Va{k) and T)ij{k) in equa tions (63) and (64) are determined by the triangle inequalities for 
k, ki and k 2 and the additional restriction of type (A2): 


'^aiX) ~ 1^^2) ^ I^1^2 ^ 2 ’ T ^2 ^ k] k k\ ^ k 2 ]k k 2 ^ j , 
'^b{k) = 1^ {ki 5^)2)^ (M^)^ I k 2 ^ ^; ki k 2 ^ k] k ki ^ k 2 T k k 2 ^ • 


They are sketched respectively in figures and 

Substituting ujiuu^ = kik 2 and = kk 2 we have 


— f 

^ ^k Jvaik) 

b _ 8^ 


L 


-Dhik) 


Si2kS{ujt’t2) in\n^2 - Hnt) ikik2f-^ dkidk2 , 
^Si2k5{ujlll) {nln\ + “ nlnl) {kik2Y~‘^ dkidk2 . 

Wi 


(67) 

( 68 ) 


(69) 

(70) 
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FIG. 6. Integration domain Vb{k) for k < 1/a/ 2 (left-hand panel) and k > 1/v^ (right-hand panel). 


Substitution of Si2k for d = 3 leads to an even simpler system, 

f dkidk2 , 

ki^k J-D^ik) 

'^k = ^ f «^2 + <«fe - nlnl) dkidk2 ■ 

^ JVb(k) 


(71) 

(72) 


E. Absence of Kolmororov-Zakharov spectra 

The central feature of wave turbulence theory is that one can compute power-law spectra satisfying the kinetic 
equations. Such Kolmororov-Zakharov (KZ) spectra are obtained for scale-invariant systems, i.e., when the dispersion 
relations and interaction coefficients are homogeneous functions of their arguments. In the KGMH model, this is the 
case only in the relativistic limit, /c ^ 1. Let us focus on this limit and exclude 2D systems: wave turbulence kinetic 
equations are known to be ill-posed in 2D for non-dispersive waves, uj ^ k^ see e.g. [12]. In 3D, the relativistic limit 
of the isotropic kinetic equations ([tT]) and ( 72 ) takes the simple form 


2'7r r 

n“ = — / {n’lnl - 2n\nl) 5{k - fci - fca) dkidk 2 , 

ki-k — J (^ 1^2 + ~ ^ 2 ^ 1 ) '^(^1 — k — ^ 2 ) dkidk2, 


(73) 

(74) 


In appendix we consider power law spectra and show that there exists no KZ solutions with local interactions 
for these relativistic kinetic equations. The only power-law solutions to these equations are the relativistic limits of 
the thermodynamic spectra corresponding to equipartition of energy (60) or of the particle invariant (61). 


We are thus left with two routes to investigate the spectral evolution of the KGMH system: in the next section we 
characterize the thermodynamic states reached in the long-time limit and study the wave condensation phenomenon. 
In sections [Vll| and [VIII| we come back to the nonlocal interactions that govern the out-of-equilibrium dynamics and we 
derive rigorously a nonlocal KZ state characterized by a nonlocal transfer of particles to a strong condensate, together 
with a local energy cascade in /c-space. This “dynamical condensation” describes the approach to the thermodynamic 
condensed state, and therefore makes the link between the two routes. 


VI. BOSE-EINSTEIN CONDENSATION 

We now study equilibration in truncated conservative systems. Non-dissipative wave systems are known to exhibit 
a condensation phenomenon when (i) there are two conserved quantities, (ii) there is a truncation at some finite 
wavenumber A, and (iii) the spatial dimension d is large enough. In discrete systems such as the array of pendulums 
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FIG. 7. The chemical potential /x as a function of the “energy per particle” T — E jM. The usual Rayleigh-Jeans spectrum is 
a solution for intermediate values of T only. For low T condensation occurs at low /c, whereas for large T energy and particles 
accumulate at the cutoff wavenumber K. 


described in section [TTj a truncation in /c-space naturally arises. In quantum systems, the Bose statistics ensures that 
modes with large enough energy remain unoccupied, hence an effective truncation. 

Condensation of classical waves was originally reported for the Gross-Pitaevskii equation mils]. Further studies 
on condensation also considered equations of the Klein-Gordon type [14]. For all these systems one can compute KZ 
solutions corresponding to an inverse cascade of the particle invariant. Condensation then naturally occurs, as the 
inverse cascade induces an accumulation of particles at low k. In contrast with these systems, we did not find any 
such KZ inverse cascade for the MHKG equation. Furthermore, the particle invariant is only an adiabatic one. In 
the following, we follow an approach similar to Connaughton et al. [7] to show theoretically and numerically that 
condensation nevertheless occurs. 

We consider decaying solutions with fixed initial number of particles JV and several values of the initial energy £. 
After some transient, we expect the system to reach the thermodynamic spectra ( [6^ , with a pair (T, /i) corresponding 
to a pair (f, A/"). However, we will see that there are some values of (f, A/") for which no pair (T, /i) can be found with 
positive fi. 

From the thermodynamic spectra, we compute the numbers Na and of a and b particles, as well as the energy Ea 
and El) contained in these two fields. Assuming for simplicity that k G [0, K] is a continuous variable and considering 
the 2D problem, we obtain 


Na = 27r L^T 
= 2 'kL‘^T 

Ni) = 2 'kL‘^T 
Ea = 2 ttL‘^T 
= 2tiL‘^T 


f 


2)i+^/2 + 


dA: 


\/2 + - V2 + 2^\ogi ^— 

^ I ^/2 + W + 2/i 


f 


0 M + 

K 


dfc = 2t^L^T 
-nr k 


2I1 + V2 + k? 


K + fiXog 


dk 






I ^ 

f + 2M^/2 - V^) + 4/l„g 


(75) 


(76) 

(77) 
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Ef, = 27rL^T 


f 


e 


k 

We define the average energy per particle as 


dk = 27 tL^T 


y (if - 2ll) + I? log (^1 + ^ 




(78) 


(79) 


T is independent of temperature, and can be specified at the beginning of a numerical simulation by an appropriate 
choice of the initial condition. However, not all values of T are compatible with the thermodynamic spectra. Indeed, 
we show in Figure a plot of /i as a function of T. There is a minimum value of T under which fi becomes negative: 
the thermodynamic spectrum is not a valid solution anymore, and wave condensation occurs. In the limit of large 
the threshold is Ter = K/3. 

We stress the fact that in the present 2D system, condensation occurs even in the infinite domain limit. This does 
not contradict the Mermin-Wagner theorem, because the dispersion relation of the 6-waves goes to zero slow enough 
k diS k ^ 0) for the integrals (75) and (76) to converge at /c = 0, even when /i = 0: for a given value of the 


energy, only a finite number of particles can be in the Rayleigh-Jeans spectrum, and the rest of the particles have 
to condense. By contrast, in the 2D Gross-Pitaevskii equation the dispersion relation is ujk = k‘^^ and the number of 
particles contained in the Rayleigh-Jeans distribution involves a diverging integral for /i = 0: for a given value of the 
energy, an arbitrary large number of particles can be in the Rayleigh-Jeans spectrum, and condensation does not occur 
in 2D. For a system with dispersion relation ^ k^ in dimension d, the criterion for Bose-Einstein condensation to 
occur is <j < d m- 

There is also a limiting value of d^(/i) that is attained for /i ^ oo. In the limit of large dC, this limit is d^(/i ^ 
oc) = K/2. As one approaches this limiting value of T, the spectrum tends to the equipartition of A7. Obviously, for 
a finite number of particles, taking the limit ja ^ oo also requires T ^ oc. If the initial values of J\f and £ are such 
that T > T{jii ^ oc), then there is an accumulation of energy in the high k modes. This accumulation corresponds 
to Rayleigh-Jeans thermodynamic spectra ( [^ with negative values of T and /i. 

To investigate numerically these different thermodynamic states, we perform numerical simulations with an initial 
condition 0k ^ and the same initial value Af — 44. Large negative values of a correspond to low T and steep 
decreasing spectra, whereas positive a correspond to large values of T. In the following we describe the stationary state 
attained after a long integration time. The size of the domain is [0,127r]^ and the resolution is 128^. After dealiasing, 
we keep 64 modes in each direction, i.e., we keep wave numbers k with k < K = 32 j6 where the step is Sk = 1/6 
for the two components of k. We plot the spectra of a and h modes, which we evaluate using |akp — cj^^^jAkP and 

|6kp — More specifically, we show the ID spectra Sa{k) and Si){k) defined in such a way that J Sa-i){k)dk 

is the variance of the corresponding field. 

The four runs described in figures to 11 illustrate respectively (i) a Rayleigh-Jeans state, (ii) the threshold for 
wave condensation, (iii) wave condensation, and (iv) accumulation at large k. In the latter simulation the system has 
not settled to a steady state: the ratio T = £/J\f very slowly decreases in time until it reaches its critical value and 
a Rayleigh-Jeans state. The drift of T probably originates from the fact that J\f is only an adiabatic invariant: it is 
not conserved by the higher order nonlinearities that we discarded in the present weakly nonlinear study. 


A. Non-uniform condensates 


To study the condensation transition, we performed a set of numerical simulations with decreasing initial values 
of T: 8iS T decreases, so does until the denominator of the Rayleigh-Jeans spectrum (62) vanishes for some value 
of k. This happens first for a 6-mode, since uj\. < uj^ for a given k. As a consequence, close to threshold we expect 
condensation of 6-waves, but no condensation of a-waves: particles accumulate in the lowest-/^ 6-modes, which in the 
equivalent particle physics system corresponds to Bose-Einstein condensation of pious. 

These 6-modes cannot have vanishing wavenumber, because 6k=o is fixed in the leading order. Indeed, inserting 
the weakly nonlinear expansion '0(x, t) = 1 + e(A(x, t) + ix(x, t)) in the expression for the vanishing charge Q, we 
obtain to lowest order that ^ (f ydx) = 0, i.e., the uniform phase mode remains constant in time, and it can be set 
to zero. As a consequence, in a finite periodic domain of size L in d dimensions, there are d 6-modes corresponding 
to the lowest non-vanishing frequency = k = 27 t/L. Eor such a finite-size system, condensation occurs as soon as 
/i = —2'k/L^ and it manifests as an accumulation of particles in these large-scale (but non-uniform) 6-modes. 

As the condensate fraction increases, the weakly-nonlinear assumption rapidly breaks down in the \ow-k modes. 
Indeed, the level of nonlinearity is determined by the ratio of the typical intensity to the typical frequency of the 
waves. Hence nonlinearity increases during condensation, because 6-particles accumulate in the low-frequency modes. 
Note that conservation of the particle invariant JV may even fail because of this increased nonlinearity. Nevertheless, 
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FIG. 8. a = —0.5: in the stationary state, T — 2.42 G (iF/3,iF/2), so we expect a Rayleigh-Jeans distribution. Top: spectra 
of a and b modes (left), and spectra normalized by a Rayleigh-Jeans spectra with fi — 1.67 and T — 1 x 10“^ (right). Bottom: 
snapshots at hnal time of the A and y helds. 





FIG. 9. a = — 0.8: in the stationary state, T — 2.12, which is close to the critical value (approximately i^/3) for wave 
condensation on this K = 32j6 grid. Top: spectra of a and b waves (left), and spectra normalized by pure energy equipartition 
spectra with T = 3.5 x 10~^ (right). Bottom: snapshots at hnal time of the A and y helds. 
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FIG. 10. a = —1.3: in the stationary state, T — 1.12, which is below the value for wave condensation. Top: spectra of a and b 
waves (left), and spectra normalized by energy equipartition spectra with T = 3.3 x 10“^. The flat part at high k is in thermal 
equilibrium, and the bump at low k indicates condensation. Bottom: snapshots at final time of the A and y fields. 



FIG. 11. a = — 2: after 6000 time units, T — 3.08, which is above the critical value for high-/c accumulation. Top-left: 
spectra of a and b waves. Top-center: spectra normalized by equipartition of A/", i.e., Rayleigh-Jeans spectra in the limit 
(T, /i) ^ (oo,oo), with T!jjL — 1.5 x 10“^. The flat part at low k is in equipartition of Af, and the bump at high k indicates 
accumulation. Top-right: spectra normalized by thermodynamic spectra with negative temperature T = — 1.3 x 10“^ and 
chemical potential /x = —10. Bottom: snapshots after integration over 7200 time units of the A and y helds. 
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FIG. 12. Condensate fraction 
b particles (□), and of 2Na + 


as a function of the energy per particle SjAf. We sho w th e condensate fraction of a particles 
Nb (A). The dashed-line is the theoretical prediction (80). 


we observe in the numerics that condensation is a robust phenomenon, that occurs even when the particle invariant 
is only an adiabatic one. As condensation proceeds, the enhanced nonlinearity arrests further accumulation and 
sharpening of the spectrum in the lowest (non-zero) frequency modes, resulting in a condensate that spreads over a 
few low-frequency modes. Because dynamical separation of the normal modes of the linear system is lost for strong 
nonlinearity, condensation of the 6-modes “spills” into the \ow-k a-modes, even though the weakly nonlinear theory 
predicts no a-mode condensation. The resulting condensation behavior is displayed in figure [Tq| 


B. Condensed fraction 


To extract the condensate fraction from the numerical simulations, we first determine the temperature of the system 
using a fit to the high-/c part of the spectra. These spectra consist of a fiat isothermal level, plus a bump at low-k 
corresponding to the condensate. By subtracting the isothermal part of the spectra before integrating over k we access 
the number of condensed particles. The condensed fraction is the ratio of the number of condensed particles to the 
total number of particles. We plot it in figure [T^ From the numerics, the threshold for condensation is estimated to be 
T ^ 2.1. We can predict this value accurately if we take into account the discreteness of the lattice of wave numbers. 
Replacing the integrals in equations ( [7^ to ( [7^ by a sum over the discrete lattice of wave numbers, we obtain a 
critical value for condensation Ter — 2.06, hence a very good agreement between the theory and the numerics. 

The condensed fraction can be predicted theoretically as follows: under the threshold for condensation, the chemical 
potential vanishes and the thermodynamic spectrum of b particles diverges for small k: slightly below the onset Tcr^ 
a macroscopic number of b particles condenses, while the a particles and the remaining b particles are 

(approximately) in energy equipartition, i.e., they have 0 thermodynamic spectra. The non-condensed particles 
carry (almost) all the energy E of the system, and their number is 2JVa + A/^ — . The ratio of these two 

quantities is Tcr^ from which we deduce the condensed fraction of b particles. 


VA") -FC.-.F 
Af Ter 


This prediction is in good agreement with the numerical simulations (see figure 12). 


( 80 ) 
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a. b. 

a{k) ^{k) 




FIG. 13. Reduced growthrate of 5-waves as a function of wavenumber k, for a-waves at wavenumber k. The dashed curve is 
(\/2 + — /c) + k~^. Left-hand panel: for 1, 5-waves grow at wavenumber k ~ Right-hand-panel: for ^ 1, 

the growthrate has two strong maxima, one Sii k 1 /k, and one at /c ~ 


VII. STATISTICAL DESCRIPTION OF THE DECAY INSTABILITY 

Kinetic equations allow to compute the out-of-equilibrium evolution of some initial distribution of waves. In 
particle physics, such an initial distribution can be given in some cosmological model, or it can correspond to particles 
appearing at t = 0 in an experiment. As we mentioned earlier, for the KGMH system the transfers governed by the 
kinetic equations can be very nonlocal in /c-space. As an example of such nonlocal transfers, we come back to the 
decay instability described in section El which we now study for isotropic ensembles of waves. 


A. Decay instability for an isotropic ensemble of a-waves 


Consider a 3D isotropic ensemble of a-waves with wavenumber k. These waves spontaneously induce some 6-waves 
through the decay instability. A key question is to know the distribution of 6-waves that is produced by this instability: 
what is the wavenumber (or energy) of the 6-waves? In the particle physics analogous system, we use the kinetic 
equation to predict the distribution of pions that arises from the decay of a-mesons with wavenumber k. 

The initial distribution is 


Na 


ni = 




6{k — K,). 


< = 0. 


(81) 

(82) 


As already mentioned in section 


VC 


this distribution is a steady solution to the kinetic equations because = 0. 
However, it is an unstable solution, as we can see by considering an infinitesimal perturbation n^(t) ^ 1. Inserting 


(92) into kinetic equation (72) and retaining only the linear terms in the infinitesimal perturbation nL we obtain 


nl = 


Na 




(ni 


■ n 


\/2+k^—/c 


) iffcex = 


V2 + — K ^/2 + K? + , 


(83) 


hl=Q 


otherwise 


(84) 
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where the condition in (83) arises from the integration domain Vi)(k). 

We write the same equation at wavenumber \/ 2 +^ — k to obtain the following 2 x 2 system of linear ODEs, 


= 


Na 


Nn. 


riy 


■ n 


n 


^/2-\-k^ — k 


— k)‘^K,\/2 + 


■\/2-\-k^ —k 

( 


)' 


ni 


■ n 


\/2+/^^ —/c 


). 


(85) 

( 86 ) 


if /c G X, and fi^ = h^ 2 ^^ 2 _/g = 0 otherwise. 

Consider an eigenmode (n^,n^ 2 ^^ 2 _fe) = where a is the growth rate. Substitution into the 

system of equations yields a 2 x 2 system in and ^^ 2 +av^ -/c’ nontrivial solutions when the determinant 

vanishes. This gives the following expression for the reduced growth rate a = aK.y/2 + k? L^/Na^ 

a{k) = -- 2+IT iffceX, (87) 

(V2T7?-fc)' 

d{k) = 0 otherwise. ( 88 ) 

Plots of d(/c) are provided in figure in the limiting cases k, 1 and k, ^ 1. For k, X is a narrow interval 
centered around and 6 -waves appear at this wavenumber: there is a strong nonlocal transfer of massive 

particles with hz 1 into massless particles with k For ^ 1, the growth rate has a maximum at k k. and 

a maximum at k = \/2 ^ n? — ^ ^ 1 : massless particles will appear preferentially at these two wave numbers, 

hence a strong nonlocal transfer from a-particles with k = n into 6 -particles with k = ^ <^1. 

These predictions for isotropic initial distributions of a-waves are qualitatively compatible with the numerical results 
of section [IV| on the decay of a single a-mode. If the initial a-mode has wavenumber k ^ 6 -waves appear at both 

wave numbers k 1 and k k (see figure]^, whereas if the initial a-mode has wavenumber k then 6 -waves 

appear at wavenumber k ^ (see figured. 


B. Reaction yield of the decay instability 


Starting from an isotropic distribution of a-particles only with wavenumber together with weak background noise 
on the a and b fields, the decay instability sets in and 6 -waves appear through nonlocal transfers. The intensity of the 
a-waves decreases accordingly, to ensure conservation of the energy and particle invariants. When 6 -waves become 
strong enough, the quadratic term in of kinetic equation ( plj ) cannot be neglected anymore. This term modifies 
the spectral distribution of a-waves, the system enters the nonlinear regime and eventually it reaches thermodynamic 
equilibrium, with both spectra of a- and 6 -waves spreading over the entire range of wave numbers. 

We can thus determine the saturated state of the decay instability using the thermodynamics described in section 


VI Indeed, the initial distribution (92) corresponds to an energy per particle 






(89) 


An a-wave can be seen as a molecule made of two 6 atoms, and the decay instability corresponds to the reaction 
a ^ 6 + 6 . In the saturated state, one can compute the yield r] of this reaction using thermodynamics: out of the 
initial JV number of 6 atoms (two per a-wave, and one per 6 -wave), what fraction 77 decays into free 6 particles? The 
reaction yield is given by the equilibrium value of and it depends on the initial wavenumber k of the a-waves. 

We now compute it in 2 D, as a function of k. 

For T > Ten 110 condensation occurs, and the yield is simply given by 


T] = 


Nbifi) 


2A6x(/i) + W(/^) 


K 


(t%) 


K jii in 


{n+k) 




(90) 


where the chemical potential fi is related to k through the relation (89), i.e., k = — 2 where T{fi) is given 

by (P. 
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FIG. 14. Reaction yield 77 for d = 2 as a function of the initial wavenumber n of the a-waves, for a maximum allowed wavenumber 
K ^ 1 (here K — 100). In the decreasing part, some fraction of the produced h particles condenses. In the increasing part 
such condensation does not occur. 


For F < Tcr^ condensation of 6-waves occur and the spectra consist of thermodynamic spectra with /i = 0, together 
with 6-particles condensed at low k. The reaction yield becomes 

_ Nhin = 0 ) 


= T 


Nb{ii = 0 ) 


— 0 ) + — 0 ) 


Ter 




(91) 


In figure 14, we plot the yield given by expressions ( [M| ) and (91) as a function of the initial wavenumber k of the 
a-particles, for a maximum wavenumber K ^ 1. In this limit, Ter — i^/3, and we observe the following behavior: the 
reaction yield is a non monotonous function of k. It is close to unity for / 1 : ^ 1, i.e., the reaction is almost complete, 
with most of the a-particles eventually decaying into b particles. In this regime a large fraction of the produced 
6-particles form a con densate. As k increases the yield decreases, with a minimum attained at the threshold of wave 
condensation Kcr = >/ ^T^r ~ ^ — 2*Fcr — 2iF/3. For k > Kcr the yield increases with k. For k = K the yield reaches 
0.5, i.e., half of the 6 “atoms” only are free 6-waves, the other half remaining in the form of a-waves. 


VIII. DYNAMICAL CONDENSATION AS A NONLOCAL KZ CASCADE 

In this section we describe the dynamical process by which condensation is achieved. We first show that ensembles 
of 6-waves at low-enough wave numbers are steady solutions to the kinetic equations. We refer to such solutions as 
“condensates”. We consider an infinitesimal perturbation of the a- and 6-wave spectra and study the linear stability of 
a condensate at k = 0. This analysis indicates that there is a rapid transfer of particles from the perturbation to the 
condensate. However, this transfer is only partial, and some marginally stable perturbations remain. We study the 
evolution of these remaining perturbations through a weakly nonlinear expansion and derive a kinetic equation that 
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governs the nonlinear evolution of the perturbation to the condensate. This kinetic equation admits a KZ solution, 
with a net energy flux towards high wave numbers. In a infinite system, this solution corresponds to a condensate 
that increases in amplitude, while the energy of the remaining waves cascades to larger and larger wave numbers. 
In a finite system, the cut-off in wavenumber arrests this cascade and the system eventually reaches thermodynamic 
equilibrium. 


A. Steady solutions with 6-waves only 


Consider the kinetic equations for = 0, and 7 ^ 0, and further restrict attention to distributions of 6 -waves 
at low-wavenumbers: we assume that there is a maximum wavenumber 1C that supports the 6 -waves, i.e., ^ 0 , 

and 7 ^ = 0 for 


and (72) 


k > JC. The evolution of this initial distribution of waves is governed by the kinetic equations ( plj ) 
At t = 0+, only the quadratic terms in may be nonzero. In the right-hand-side of equation (71), the 


quadratic term involves triads of wave numbers (/ci, /c2, k = ) such that 7 ^ 0 , 7 ^ 0 , kik 2 > 1 / 2 , 

ki k 2 > \/2. Such ki and k 2 cannot be found if /C < ^: provided the distribution of 6 -waves is restricted to wave 
numbers smaller than the right-hand-side of 

0 vanishes. We assume that this condition is satisfied by the 


initial distribution, and we consider the right-hand-side of kinetic equation (72). Because = 0, only the quadratic 
term in may be nonzero. However, it is nonzero only if we can find two wave numbers k and k 2 such that 7 ^ 0, 
7 ^ 0 , and kk 2 > 1 / 2 . But the latter condition cannot be satisfied if /C < hence the right-hand-side of (72) 
vanishes. 

To summarize, we have shown that an initial distribution of 6 -waves only = 0) is a steady solution to the kinetic 

equations (71) and (72) provided these waves have only wave numbers smaller than 


B. Dynamics of waves on a condensate 


For simplicity, we consider a 6 -wave condensate at very low-wavenumber, which we simply write as A5(k)^ and we 
denote as n^{k) and n^{k) the remainders of the wave spectra. Hence, 

nl=nl, (92) 

nl = AS{k) + hi . (93) 

Inserting this decomposition into the 3D kinetic equations for /c 7 ^ 0 leads to 
jia 4:71 A /i, 271 


riu = 


kV2 + fc 2 

AttA 


J JVaik) 


2fi[h^) dkidk2 . 


(^Vk^ “ + TI / A'^(^i2k) {^>2 + nlnl - nlnl) dkidk2 . 


(94) 

(95) 


The terms proportional to A arise from the interaction between the condensate and the remainders, while the integral 
terms correspond to interactions of the remainders. The latter terms are exactly the collision integrals of the ori gina l 
kinetic equations, with n replaced by h. We stress the fact that if and are the isothermal distributions ( [ 6 Q| ), 
we obtain a steady solution to (94) and ([^, i.e., both the terms proportional to A and the integral terms vanish: a 


condensate at k = 0 together with remainders in isothermal equilibrium constitute a steady solution to the kinetic 
equations, in perfect agreement with the thermodynamic treatment of Bose-Einstein condensation. 


To study the dynamical process of condensation, we perform a multiple-scale expansion of equations (94) and (95), 


to describe the evolution of remainders much weaker than the condensate, n A. Consider a small parameter e <C 1 
and A = 0 ( 1 ), and expand the remainders as 

(96) 

(97) 


nz = en 


“’°(fc,i,T)+ e2 n'^'d{k,t,T) 


n\ = 


{k,t,T) + e^h^'^{k,t,T) + ... , 


where the slow time is T = et. 


Inserting this decomposition into (94) and ( [9^ gives at order 0(e) the linear stability analysis for the condensate, 

47rA 


dtnf = 


O -6;0 47r74 /_a;0 


- 6;0 


If) 


(98) 




;») 


(99) 
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Considering the first of these two equations at wavenumber \/k‘^ — 2 instead of k leads to the following 2x2 system 
of first-order linear ODEs, 




AttA 


~6;0 47rA / a;0 ~6;0\ 

{^^2 -^k ) ■ 


yJkP' — 2 k 


( 100 ) 

( 101 ) 


Solutions to this system of ODEs can be sought in the form ^ ^ which leads to the following growth 

rates, 




< 0 . 


Decomposing on the corresponding eigenmodes we obtain 

h^’° = Ci{k,T) - . 

Because (72 < 0, the contribution in C 2 rapidly decays on the fast time t, and after a short transient, 


“ a;0 


= hf = nfe(T), 


( 102 ) 

(103) 

(104) 

(105) 


One can check that the C 2 part of the spectra has vanishing total energy at any time. Its rapid decay occurs 
at constant energy, as it should. By contrast, this C 2 contribution carries a nonzero number of particles A /2 = 
/ C 2 {k)y/k‘^ — dk. The decay of the C 2 contribution on the fast time t therefore corresponds to a rapid 

transfer of particles between the remainders and the condensate. These rapid transfers of particles ensure that (105) 
is satisfied adiabatically. The dynamics of the system is then characterized by a single slowly-evolving spectrum 
n/c(T) from which both spectra of a and b particles can be deduced. We pursue the expansion to determine the kinetic 
equation that governs the evolution of nk{T). 

To order O(e^), we consider a linear combination of equation (94), considered at wavenumber — 2, and of 
equation ([9^, to obtain 


dt 


k‘^ -2 + k ) = - (fc + - 2^ drnk 


(106) 


27r 

+ T 


'^ki'^k — ki 77 //- dki 


[ ^ ^ 

/fc-Vfc^-2^7 ^fc + Vfc^-2 

J -- 

I '^k\'^k\—k T '^ki—k'^k dki 

J /ci >/c+ -A 


To obtain this form for the second integral, we have made the change of variables y^f + 2 ^ ki in the integrand. 
The solvability condition demands that the right-hand side vanish, which leads to the following kinetic equation. 


dr'^k — 


27r 




k-Vk^-2 


<ki<- 


k (^k \/k‘^ — 2) 

T 2 / HkiTlki—k T '^ki'^k '^ki—k^k dk\ 

Jkr>k+^ 


^ki^k — ki 277//j;^ TTz/j 


(107) 


i//v-r 2k 

In the relativistic limit k ^ 1 the equation becomes scale-invariant 




j ^ki^k — ki ^TTz/j^77//- dk\ T 2 / T Tlki^lk "^ki—k^k dk\ 

Jki<k Jki>k 


(108) 


This equation conserves the total energy, but it does not conserve the number of particles in the remainders. Indeed, 
these remainders rapidly exchange particles with the condensate, and only the total number of particles contained in 
both the remainders and the condensate is conserved. 
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In appendix I C 2| we determine the steady power-law solutions to this equation ^ As already mentioned, one 
such solution to the kinetic equation is the isothermal spectrum rik ^ k~^ that coexists with the condensate when 
wave-condensation occurs. But this kinetic equation also admits a solution of the KZ type, with exponent x = — 
This turbulent spectrum carries a flux V of energy through the scales of the system. It reads (see appendix C 2 for 
details) 


rik = 


V2k-i 


47r\/27r — 8 In 2 


(109) 


In this cascading state, the flux of energy is due to local interactions in /c-space, while nonlocal interactions continuously 
transfer particles from the remainders to the condensate to satisfy the adiabatic constraint rik = = ri\ (one could 

highlight the corresponding slow growth of the condensate by considering the 0{e^) equations at /c = 0 and introducing 
an evolution of the amplitude A of the condensate on an even slower timescale e^t). 


IX. DISCUSSION 


We have developed a weak turbulence theory for the Klein-Gordon equation with a Mexican-hat potential 0 , a 
universal nonlinear equation which displays spontaneous symmetry-breaking. The fundamental state of this equation 
is a uniform field at the minimum of the Mexican-hat potential. The weak-turbulence theory describes the interaction 
of random weakly nonlinear waves propagating over this uniform background state. Such wave disturbances are of 
two-kinds, with two branches in the dispersion relation: in the particle physics context, these branches correspond to 
a massive and a massless particle fields, and the KGMH equation provides the simplest particle physics model for a 
system of interacting sigma mesons and pions. 

In the weak-nonlinearity limit, the dominant interaction between these fields is a three-wave process by which two 
massless waves can merge into a massive one, or inversely one massive wave can decay into two massless ones. This 
process conserves energy and an equivalent of the number of particles. 

We have derived the wave turbulence kinetic equations describing the coupled evolution of the spectra of massive 
and massless waves under this three-wave interaction. They admit a two-parameter family of thermodynamic equi¬ 
librium spectra of the Rayleigh-Jeans type. When the initial condition corresponds to little energy per particle, the 
thermodynamic state reached in the long-time limit displays massless-wave condensation, the classical analogue of 
Bose-Einstein condensation. 

The kinetic equations admit no Kolmogorov-Zakharov (KZ) cascading states in the usual sense. Instead, the out- 
of-equilibrium dynamics is dominated by nonlocal interactions in scale-space. We have highlighted two cases of such 
nonlocal evolution. The first one is when a massive mode spontaneously decays into massless waves, and the second 
one is when small-scale a and b fields evolve in presence of an intense large-scale b condensate. In the latter regime 
some particles are rapidly transferred to the condensate, and the small-scale remaining waves follow a single reduced 
kinetic equation which admits a KZ solution corresponding to an energy cascade. 

These different regimes are illustrated in the 3D numerical simulation reported in figure!^ the domain is [0,1/]^ 
with!/ = 27r \/6 15.4, with a resolution of 192^. The wave vector k is on the 3D grid [—96,96]^, but after 

dealiasing we keep only the modes with k < ^ 20. We plot the ID wave spectra Sx{k) = (^)^ 47r/c^|A/cp and 


'S'xW = (^) 47rA:2|xfep. 

The simulation starts off with a distribution of massive waves only, together with a weak level of massless waves. 
The initial spectrum Sx is flat up to wavenumber /c = 6, i.e., the massive waves are mostly relativistic (see panel [T^). 
The initial evolution of this spectrum corresponds to a fast decay instability. As discussed in section 1^ the decay 
of relativistic massive particles with wavenumber ka produces massless waves with ki, ka and with ki, 1. The 
latter waves accumulate in a strong condensate of massless particles. The subsequent evolution corresponds to waves 
superposed to a strong condensate, the dynamics of which we presented in section [VIII| the spectra of the small-scale 
remaining a and 5-waves rapidly become identical (see panel T^), and follow the reduced kinetic equation (108). 
After some transien t, we obtain the quasi-static spectral distribution plotted in panel [T^, which corresponds to the 
KZ energy cascade ( |1Q9[ ). Indeed, we computed the energy leaving the domain k e [1; 10] per unit time to determine 
the v alue of the energy flux. We obtained V — 7. 10“^. In the relativistic regime, the KZ spectra corresponding to 
(109) are given by kSx{k) = kS^{k) = 47r/c^n/c, i.e.. 


Sx{k)=S^{k) = 


\/27r — 8 ln2 


( 110 ) 


This theoretical prediction is plotted in figure [T^. It has no fitting parameters and is in good quantitative agreement 
with the spectra obtained numerically in the inertial range. This confirms that the system displays a nonlocal 
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a. 



b. 



c. 



FIG. 15. Out-of-equilibrium spectra Sx (blue solid line) and Sx (red dashed line) in a 3D numerical simulation of the KGMH 
equation. Panel a shows the initial condition {Sx{k) ~ 2 10“^ is out of the figure). Panel b corresponds to t = 40: after a fast 
transient, part of the massive waves decayed into massless waves. A strong condensate of 5-waves appear, and the small-scale 
remaining spectra follow the adiabatic constraint Sx — Sx- Panel c cor responds to t = 475: energy cascades from low to high 
wave numbers, in quantitative agreement with the KZ prediction (110), shown as a thick solid line. 
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Kolmogorov-Zakharov cascading state, where nonlocal transfers of particles with the condensate coexist with a local 
energy cascade. This state is therefore the intermediate asymptotics that describes the evolution between the initial 
distribution of waves and the thermodynamic equilibrium attained in the long-time limit. 
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Appendix A: Three-wave resonances 


The interaction Hamiltonian (26) corresponds to three-wave processes of type a + a + a 
b 


0, CL CL —y CL, CL b b 


0,a b b and b b a. Out of these processes, the only ones allowed are those that can satisfy both the 
wavenumber and the frequency resonance conditions. 
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1. a + a + a^O and a b b ^ 0 processes 

Clearly, these processes cannot be resonant because the sum of three positive frequencies cannot be zero. Thus we 
neglect the corresponding terms in the Hamiltonian. 


2. a-\- a ^ a process 

The frequency resonance condition for this process is After the substitution ks = ki + k 2 , it 

becomes 


'\J2 k\ + y^2 + — a/ 2 + (ki + k2)^. 

Taking the square of both sides, we obtain 

k\ ‘l'\J (2 + ^i)(2 + k‘ 2 ^ = 2 k\ k^ 2(ki • k2). 


or 


1 + y^(2 + k\){2 + — (ki • k2). 

This equation has no solutions because the left hand side is strictly greater than kik 2 whereas the right hand side is 
less than kik 2 . Thus we neglect the corresponding terms in the Hamiltonian. 


3. a + b ^ b process 

For this process the resonance conditions give 

^2 + kl +k2 = y(ki+k2)2. 

Squaring both sides yields 

2 ki k2 2k2^2 k‘^ = k\-\- k^-\- 2(ki • k2). 


or 

1 + k2'\J2 k\ = (ki • k2). 

Again, this equation has no solution because the left hand side is strictly greater than kik 2 whereas the right hand 
side is less than kik 2 ^ so we neglect the corresponding terms in the Hamiltonian. 


4. b -\-b ^ a process 


For this process the resonance conditions give 

ki^k2 = a/ 2 + (ki + k2)^. 


Squaring both sides, we obtain 


k\ k^ 2k2ki — 2 ki k2 2(ki • k2), 


or 


2 k 2 ki = 2 + 2(ki • k 2 ) = 2 + 2 k 2 ki cos 0, 
where 0 is the angle between ki and k 2 . We have 


— 1 < cos 0 = 


k2ki - 1 


< 1, 


(Al) 
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i.e., such three-wave resonances exist if 


hh > 


(A2) 


This condition predicts nonlocal interactions if the initial spectrum of a-waves is at /c ^ 1: then the initial evo lution 
generates b-waves at k ^ Indeed, k ^ 1 implies ki ~ k 2 and cos^ ~ —1, which according to (Al) gives 

ki ^ k 2 ^ IIV 2 . Conversely, two short b-waves can produce a long a-wav e thr ough a wave-beating process. 

On the other hand, if the a-particles are relativistic then the condition (A2) is not restrictive: for k,ki,k 2 ^ 1 we 
have cosO ^ 1, which is typical of acoustic-like systems: the wave vectors are almost colinear in all resonant triads. 


Appendix B: Wave turbulence derivation 

Wave turbulence deals with statistical ensembles of waves. For such wave fields, a small amount of dispersion 
rapidly leads to the complex amplitudes of the different waves being uncorrelated. We can therefore consider the 
amplitudes and the phases of the waves to be random independent variables to start with. Such a random phase 
and amplitude approximation is a useful shortcut to obtain the wave kinetic equations. The brackets (...) denote an 
ensemble average over these random phases and amplitudes. 


1. Separating the time scales 

Consider the interaction representation variables 

ctk = Ck, &k = d^, 

where the small parameter e <C 1 has been introduced to study the weakly-nonlinear behavior. These variables are 
steady in the absence of nonlinearity, and they evolve slowly for weak nonlinearity, with 


Ck = -ie y] V^2 did 2 e , 


ki,k2 


dk = -2ie y] y2kCid2e®‘‘'“2^5ij^^ 

ki,k2 


(Bl) 

(B2) 


where 


abb _ , ,a b 


2 : 


abb _ , ,a b b 
‘^12k — “ ^2 “ 


Consider a time T which is intermediate between the (short) period of the linear oscillations and the (long) nonlinear 
time-scale, 

27r 27r 

uj e^uj 

where the inequality must be satisfied for both uj = and uj = uj^. 

Seek the solutions at time T via a regular expansion in e 

Ck(T) = -h -h + ..., (B3) 

dk{T) = 4°^ + + • • • • (B4) 

Below we neglect the orders of e higher than two. Substitution into the evolution equations (Bl) and (B2) gives at 
leading order = 0, i.e.. 


Ck ^ = Ck(0), d]^’ = dk(0). 


7(0) 


At the next order we have 


4'^ = -* E 

ki,k2 

4 ) = - 2 i ^ 

ki,k2 


(B5) 

(B6) 

(B7) 
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where 


= / e* 
Jo 


T ixt _ 1 

^- 

0 


Doing one more step of the recursion, we get 

4“’ = -2> E 

ki,k2 

4” = -2. E < i4‘>4“>- + 4”'4“-)e‘“S‘«k. 


^p-^^ki2^Ak 
^ ^125 


ki,k2 


Substituting here and from |B 6 | and (B7) and integrating over time yields 


^(2) _ 4 T/kl/S ^(0)^(0)^(0)* ,a66 ,,a66Nr3 rk 

Ck ——4 2_^ Vl2^A2di C 3 (24 4l/(CJ342,CJkl2j^42^12: 

ki,k2,k3,k4 

42 ) = 2 ^ y 2 i [-14 df -C^?2t)44 


ki,k2,k3,k4 


where 


+2Vi -a;?2'^)4]4k- 

i;(a:,y)= At{x)e-^y* dt, 

Jo 


<^341 = CJ3 — CJ4 — CJi, etc. 

Introduce the wave intensities 

4 = |ckp = e^l^kp, 4 = |4p = 

and expand them to second order in e, 

Iii{T) = |4°^ + = |4°y + e(4°^*4^^ + c.c.) 

+e2|4^y + e^(Ck + c-c.), 

Jk(T) = |4°^ + = Mk^l^ + + c-c.) 

+e^l4y + ^ + C-C-) • 


(B8) 

(B9) 

(BIO) 

(Bll) 

(B12) 


(B13) 

(B14) 


2. Phase averaging 

We consider the phases of the waves as ra ndom variables that are uncorrelated for different wave vectors. When 
phase-averaging the expressions (B13) and (B14), the linear order in e drops out, as usual in Wave Turbulence. 
Denoting the average over those phases as (.. .) 0 , we obtain 


(l4f)0 = 2 E 

ki,k2 

ki,k 2 

where = |Ck^^p and Jk^^ = Similarly, 

= -4 E ^ i 2 v !2 <1)442, 

ki,k2,k3,k4 

= -4 E in^2l"4°4f^«ta;gf2)42, 

ki,k 2 


(B15) 

(B16) 


(B17) 
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= 2 E [-^34 -^?2k)44 

ki,k2,k3,k4 

+2V4I (c(°)cf *df df -cv^^^)sl2]sk 

= 4 ^ \Vl\^ [/f^ - 4°^]jL°^E{-wf^l -cofl)5l^. 


(B18) 


ki,k 2 


We can further average over random initial amplitudes: the intensities at different wave numbers are independent 
random variables, for instance (4^^ 4^^) = (4^^) (4^4 


3. L ^ 00 and T ^ 00 limits. Kinetic equations. 


In wave turbulence, the limit L ^ oc is taken before the limit T ^ 00 , the latter being equivalent to the e ^ 0 
limit. 

For the limit L ^ 00 we follow the standard rules to switch to a continuous description. 



(iki(ik 2 . 


and replace the Kroenecker deltas with Dirac deltas: 

Kroenecker-^ ^ 

For the T ^ 00 limit, we take into account the following limits. 



X Dirac-(5. 


\AT{x)f 2 ttT5{x) 


and 


^[E{x^x)\ 7rT6{x). 

Introduce the wave spectra, 

"“k = (s) <''•)' = (i) ■ 

Using {Ik{T)) = (//c(0)) + T ^4^? {Jk{T)) = (7/^(0)) + T ^4^? and performing the infinite-domain and long-time 
limits, we arrive at the kinetic equations 

fik = 47r J iqip {n’inl - 2n5nk)d(wki2)d^2 dkidk2 , (B19) 

fik = Stt J |l/2kl^ {'iT'1'4 + - «2«k)<^(wi2k)<^k dkidk2 . (B20) 

Recall that here ^5^2 5 denote Dirac deltas rather than Kroenecker deltas. 


Appendix C: Isotropic power-law spectra in the relativistic limit 
1. Calculation of the collision integral for power-law spectra 


We consider the case d = 3, and we assume power-law spectra of the form 

nf = Ca,bk^ . (Cl) 

At first sight, we may be tempted to consider a more general form in which the index is different for and n^, but 
that proves unfruitful: different indices would mean different degrees of homogeneity of the different contributions 
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FIG. 16. Collision integrals as a function of the index x. a: ^^^ 2 k 2 x-i ; b: ^^^ 2 k 2 x-i ; c: ^ b: • 


to the rhs of the kinetic equations and they could not sum up to zero. Secondly, to find the steady-state power-law 
solutions, one might like to use the so-called Zakharov transformation, as is customary in WT theory EH]. However, 
we will see below that in our system all integrals can be evaluated analytically for any exponent x, and, therefore, 
the Zakharov transformation is unnecessary. 

Inserting (Cl) in equation (73) and integrating over k 2 ^ one gets 


ott r 

[Gfci (fc - kiT - 2CaCbk^k^^]dk 

Changing variable to u = ki/k, we obtain 

nl = 27Tk^^-^ Cl [ u^{l-ufdu-2CaCb [ 

L Ju =0 Jo 


1 • 


^du 


Performing the same steps to equation (74) yields 


POO POO 

47rfc2®-i CaCbJ {u^{u-lf +u=^)du-Cl J {u-lfdu 


(C2) 


(C3) 


(C4) 


The convergence of these integrals depends on the value of x, Ca and Cb- For x > 1/2 the rhs of (C4) diverges unless 
Ca = Cb- For X g]0; 1 /2] and x G [—1/2; 0[ the rhs of (C4) diverges. For re = 0 it converges only if Cb = 2Ca, for which 
it vanishes. The rhs of (03) also vanishes for (rr = 0, C b = 2Ca)- The corresponding solution is the equipartition of 
particl e nu mber M. For x g] —2 ;—1/2[ the rhs of ( |C4| ) converges provided Ca = Cb- When this condition is met the 
rhs of (C3) also converges, and the equations become 

F 2 (l+rr)' 


hi = 27rC'2)!c2*-i 


+ 


1 + re F(2 + 2rc) 


hi = 2TrC'^k 


2 7 , 2 ®-! 


F(-l-2rc)F(l 


1 X 


F(-rr) 


(C5) 

(C6) 
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where C = Ca = C^. We plot in figure the arguments of the square brackets in the right-hand side of these 

two equations, i.e., we plot respectively 27 ^ 0 ^ have a single zero at x = —1 for the range 
X g] — 2; —1/2[: this can be either the equipartition of wave energy, or the KZ spectra for a cascade of M. 

It is insightful to also plot the linear combinations 2 ^ 0 ^]^^-^ indicates if the spectrum 

of energy is stationary, whereas the second one indicates whether the spectrum of particle number Af is stationary. 
We see that fi^ + vanishes at x = — 1 but also at x = —3/2. On dimensional grounds (see e.g. [4]), one can see 
that the latter value corresponds to the KZ spectra for an energy cascade. However, this spectrum is not a stationary 
solution to the full system of equations because and do not vanish independently at x = —3/2. 

The combination 2h^ has a double zero at x = — 1. Indeed, this index corresponds both to energy equipartition 
and, on dimensional grounds, to a constant flux of Af. Still it remains unclear whether it is possible to observe a 
stationary state with such a nonzero flux of JV: the value x = — 1 is a single zero of which corresponds to 

energy equipartition. It is thus possible that is not zero when there is a flux of Af. 


2. KZ cascade on a condensate 


The remainders that evolve together with a strong condensate follow the kinetic equation (108) in the relativistic 


limit, where the a- and 6-wa ves h ave the same spectrum nk{T). Let us look for power-law solutions Uk = Ck^ to this 
equation. Substituting into (108) and changing the integration variable to u = ki/k leads to 

drUk = T:C‘^k^^~^J{x) , 


where 


nl n + 00 

J{x) = / - uf - 2 u^du ^2 / u^{u - 1)^ -{u- If du 

4 2r(-i-2x)r(i-hx) r2(i-hx) 


(C7) 


(C8) 


1 X 


n-x) 


r(2 + 2x) ' 


The function ^{x) vanishes for x = — 1 and x = —3/2, the former corresponding to energy equipartition and the 
latter being a KZ-type solution. Let us compute the energy flux V in this KZ state: the kinetic equation is a 
conservation equation for the total energy, that we can write in terms of the energy spectrum {uJa + ujfd'KkAuk = 
Sirk^Uk. Considering power-law solutions rik = Ck^, 


driSTTk^Uk) = C‘^J{x ). 

OjK 


(C9) 


The flux V of energy going from k < ko to k > ko per unit time is obtained by integrating the right-hand-side 
expression from k = 0 to ko, 

" -,J{x). (CIO) 


2^2 ^0 


V = -S7r^C' 

2x T 3 

For the KZ exponent x = —3/2, this flux is independent of /cq, indicating a con stant flux of energy cascading through 
the scales of the system. Indeed, both the numerator and the denominator of (CIO) vanish at x = —3/2, and Taylor 
expansion of ^{x) to first order around this value leads to 7^ = 327r^(7^(7r — 4 In 2). We can therefore express C in 
terms of V, which leads to the KZ solution, 

. 3 

(Cll) 


rik = 


V^k- 




Tn2 


Appendix D: Numerical code 


Let us rewrite the complex Tachyonic Klein-Gordon equation as 

V’tt - At/) +(-5 += 0 . (Dl) 

We want to solve this equation in a periodic domain [0,27r]^. The parameter g > 0 allows to put the crossover 
wavenumber anywhere in k space, i.e. the dispersion relation for a-waves is now cjg = ^/2g k‘^. For better accuracy, 

one wants to integrate exactly the rapidly oscillating part of the solution. Here we write the equation as: 

ftt = + AkL with NL = (IV^I^ — g)f . 


(D2) 
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If we drop the term NL^ a solution in Fourier space is iIj\^ = . Let us use the method of variation 

of constants: we write and impose the constraint A'j^{t)e^^^ + A'_{t)e~'^^^ = 0. Inserting 

this decomposition into equation (D 2 ) we obtain 


A'_ = -^e^^*NLy,, 


(D3) 

(D4) 


where NL\^ is the Fourier amplitude of 


g)'ip for wavevector k. We integrate equations (D3) and (D4) using an 


Adam-Bashforth scheme with time-step dt: 


1 di^ 

A+{t + dt) = A+{t) + [3NLk{t) - - dt)] , 

1 di^ 

A_(t + dt) = A_{t) + \-3NLk{t) + e-*^‘^*NLk{t - dt)] . 

Ak ^ 


Let us define and (t) = A_{t)e 


—ikt. 


{t + dt) = e 
ip^ (t + dt) = e 


ikdt 


i dt 


ikdt 


—ikdt 


^i{t) + — {3NL^it)-e 
d) + ^ {-3NLk{t) + e 


NL^{t-dt)) , 

-ikdtNL^{t - dt)) , 


(D5) 

(D6) 


(D7) 

(D8) 


For nonzero the Fourier amplitude of the field is then '0k(^) = V^k (0 + (^)' mode k = 0 must be taken care 

of separately. Using a standard finite-difference approximation to the second order time-derivative, we get 

'0o(^ + dt) = 2tpo{t) — 'ipo{t — dt) — dt^NLo . (D9) 


We compute the nonlinear term using a standard pseudo-spectral method with dealiasing following the 1/2 rule. A 
time step consists of the following: the nonlinear term is computed in real space. The field is computed as the 
inverse Fourier t rans fo rm o f a nd t he nonlinear term is transformed back to Fourier space to get NL\^{t). 

Then equations (D7), (D8) and (D9) are iterated to obtain at time t + dt. 
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